Packages / R-options:

.libPaths('~/R/x86_64-pc-linux-gnu-library/4.1/')
suppressPackageStartupMessages({
  library(plotly)
  library(gaston)
  library(dplyr)
})
# load engine's functions
sapply(FUN = source,
       X = list.files("../../src",
                      pattern = ".R$",
                      full.names = T))

#  R options
options(max.print = 20)

1 Get data

1.1 Inputs

genoFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/phasedGeno_phased.vcf.gz'
crossTableFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/crossTable.csv'
SNPcoordFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/linkageMap.csv'
markerEffectsFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/markersEffects.json'

1.2 Predictions

predictions <- calc_progenyBlupEstimation(
  genoFile = genoFile,
  crossTableFile = crossTableFile,
  SNPcoordFile = SNPcoordFile,
  markerEffectsFile = markerEffectsFile,
  # outFile = outFile
)
predictions$Cross <- paste0(predictions$ind1, '_X_', predictions$ind2)

1.3 Crossing simulation

set.seed(1234)
markerEffects <- readMarkerEffects(markerEffectsFile)
crossSimOutFile <- crossingSimulation(
  genoFile = genoFile,
  crossTableFile = crossTableFile,
  SNPcoordFile = SNPcoordFile,
  nCross = 100)
simulation <- gaston::read.vcf(crossSimOutFile)

## Warning in gaston::read.vcf(crossSimOutFile): NAs introduced by coercion

## Warning in gaston::read.vcf(crossSimOutFile): Some unknown chromosomes id's (try
## to set convert.chr = FALSE)

simulation <- gaston::as.matrix(simulation)
simulation <- simulation[, row.names(markerEffects$SNPeffects)]

simulation <- data.frame(
  simBV = as.vector(simulation %*% as.matrix(markerEffects$SNPeffects)) + markerEffects$intercept,
  Cross = as.factor(sub('-\\d+$', '', row.names(simulation)))
)

2 Simulation VS Predictions

p <- plot_ly() 

# add boxplot of simulated BV
p <- p %>% add_boxplot(data = simulation,
                       x = ~Cross,
                       y = ~simBV,
                       name = "Boxplot of Simulated individuals' BV")
# add jittered markers of simulated BV
p <- p %>% add_boxplot(data = simulation,
                       x = ~Cross,
                       y = ~simBV,
                       line = list(color = 'rgba(0,0,0,0)'),
                       marker = list(color = 'rgba(31, 119, 180, 0.5)'),
                       fillcolor = 'rgba(0,0,0,0)',
                       name = "Simulated individuals' BV",
                       boxpoints = "all",
                       pointpos = 0,
                       jitter = 1,
                       hoverinfo = 'text',
                       text = apply(simulation, 1, function(l) {
                         paste(names(l), ":", l, collapse = "\n")
                       }))


# add predicted BV
p <- p %>% add_markers(
  data = predictions,
  x = ~Cross,
  y = ~blup_exp,
  marker = list(color = 'rgb(255,0,0)'),
  name = "Theoretical BV",
  error_y = ~ list(array = 2*sqrt(blup_var),
                   type = 'data',
                   color = '#000000'),
  hoverinfo = 'text',
  text = apply(predictions, 1, function(l) {
    paste(names(l), ":", l, collapse = "\n")
  })
)

p

2.1 Test observed mean

alpha <- 0.05
beta <- 0.95
pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>% 
  filter(Cross != "Coll0659_X_Coll0425") %>% 
  group_by(Cross) %>%
  summarise(obs_mean =  mean(simBV),
            blup_exp = unique(blup_exp),
            obs_var = var(simBV),
            blup_var = unique(blup_var),
            p.value = t.test(x = simBV, mu = unique(blup_exp))$p.value,
            delta = power.t.test(n = 100, power = beta, sd = sqrt(unique(blup_var)), sig.level = alpha)$delta)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim

idLine <- data.frame(mean = c(min(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) - 1,
                                max(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) + 1),
                     var = c(min(c(pred_vs_sim$blup_var, pred_vs_sim$obs_var)) - 1,
                                max(c(pred_vs_sim$blup_var, pred_vs_sim$obs_var)) + 1))

rmse <- sqrt(mean((pred_vs_sim$blup_exp - pred_vs_sim$obs_mean)^2))

RMSE = 0.6766344

#plot mean
plot_ly(type = "scatter",
        mode = "markers",
        data = pred_vs_sim,
        x = ~blup_exp,
        y = ~obs_mean,
        name = "Observed mean vs prediction",
        hoverinfo = 'text',
        text = apply(pred_vs_sim, 1, function(l) {
          paste(names(l), ":", l, collapse = "\n")
        })
) %>%
  add_lines(inherit = FALSE,
            x = idLine$mean,
            y = idLine$mean,
            name = "Identity line")

2.2 Test normality

pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>% 
  filter(Cross != "Coll0659_X_Coll0425") %>% 
  group_by(Cross) %>%
  summarise(obs_var = var(simBV),
            blup_var = unique(blup_var),
            p.value = ks.test(simBV, "pnorm", mean = unique(blup_exp), sd=sqrt(unique(blup_var)))$p.value)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim

#plot var
plot_ly(type = "scatter",
        mode = "markers",
        data = pred_vs_sim,
        x = ~blup_var,
        y = ~obs_var,
        name = "Observed variance vs prediction",
        hoverinfo = 'text',
        text = apply(pred_vs_sim, 1, function(l) {
          paste(names(l), ":", l, collapse = "\n")
        })
) %>%
  add_lines(inherit = FALSE,
            x = idLine$var,
            y = idLine$var,
            name = "Identity line")

Appendix

Session Information (click to expand)

## Document generated in:
## Time difference of 27.34104 secs
## 
## CPU: AMD Ryzen 5 3600X 6-Core Processor
## Memory total size: 32.79236 GB
## 
## 
## Session information:
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Pop!_OS 22.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_1.0.10       gaston_1.5.9       RcppParallel_5.1.6 Rcpp_1.0.9        
## [5] plotly_4.10.1      ggplot2_3.4.0     
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0    xfun_0.36           bslib_0.4.2        
##  [4] memuse_4.2-2        purrr_1.0.1         splines_4.1.1      
##  [7] lattice_0.20-44     colorspace_2.0-3    vctrs_0.5.1        
## [10] generics_0.1.3      htmltools_0.5.4     viridisLite_0.4.1  
## [13] vcfR_1.13.0         mgcv_1.8-36         yaml_2.3.6         
## [16] utf8_1.2.2          rlang_1.0.6         jquerylib_0.1.4    
## [19] pillar_1.8.1        glue_1.6.2          withr_2.5.0        
## [22] DBI_1.1.3           lifecycle_1.0.3     stringr_1.5.0      
## [25] munsell_0.5.0       gtable_0.3.1        htmlwidgets_1.6.1  
## [28] breedSimulatR_0.3.2 evaluate_0.20       knitr_1.41         
## [31] permute_0.9-7       fastmap_1.1.0       crosstalk_1.2.0    
## [34] parallel_4.1.1      fansi_1.0.3         pinfsc50_1.2.0     
## [37] scales_1.2.1        cachem_1.0.6        vegan_2.6-4        
## [40] jsonlite_1.8.4      digest_0.6.31       stringi_1.7.12     
## [43] grid_4.1.1          cli_3.6.0           tools_4.1.1        
## [46] magrittr_2.0.3      sass_0.4.4          lazyeval_0.2.2     
## [49] tibble_3.1.8        cluster_2.1.2       ape_5.6-2          
## [52] tidyr_1.2.1         pkgconfig_2.0.3     ellipsis_0.3.2     
## [55] Matrix_1.5-3        MASS_7.3-58.1       data.table_1.14.6  
## [58] assertthat_0.2.1    rmarkdown_2.19      httr_1.4.4         
## [61] rstudioapi_0.14     R6_2.5.1            nlme_3.1-152       
## [64] compiler_4.1.1

shiny::HTML('&nbsp;
<!-- Add icon/font library -->
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.13.0/css/all.min.css">
<link rel="stylesheet" href="https://fonts.googleapis.com/css?family=Lato">

<div class="footer" style="font-family: Lato">
    <hr />
    <p style="text-align: center;"><a href="https://github.com/juliendiot42">Julien Diot</a></p>
    <p style="text-align: center;"><span style="color: #808080;"><em>juliendiot@ut-biomet.org</em></span></p>

<!-- Add font awesome icons -->
    <p style="text-align: center;">
        <a href="https://github.com/juliendiot42" class="fab fa-github"></a>
        <a href="https://www.linkedin.com/in/julien-diot-949592107/" class="fab fa-linkedin"></a>
        <a href="https://orcid.org/0000-0002-8738-2034" class="fab fa-orcid"></a>
        <a href="https://keybase.io/juliendiot" class="fab fa-keybase"></a>
    </p>
</div>
&nbsp;')

Crossing simulation vs Predictions

Julien Diot

2023/01/27

Packages / R-options:

1 Get data

1.1 Inputs

1.2 Predictions

1.3 Crossing simulation

2 Simulation VS Predictions

2.1 Test observed mean

2.2 Test normality

Appendix